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The nonlinear response of symmetric and loaded 
airfoils to an impinging vortical gust is investigated in 
the parametric space of gust dimension, intensity, and 
frequency. The study, which was designed to 
investigate the validity limits for a linear analysis, is 
implemented by applying a nonlinear high-order 
prefactored compact code and comparing results with 
linear solutions from the GUST3D frequency-domain 
solver. Both the unsteady aerodynamic and acoustic 
gust responses are examined. 

INTRODUCTION 

Most unsteady aerodynamic and aeroacoustic analyses 
considering interaction of upstream flow 
nonuniformities with downstream structural 
components follow linear theory' wherein the unsteady 
part of the flow is considered a small perturbation about 
the steady mean flowfield. This approach, which relies 
on solutions to the linearized Euler equations, is 
adopted, for example, in problems dealing with rotor- 
stator interactions where the rotor wakes are 
represented as low-intensity flow disturbances (vortical 
gusts) convected by the mean flow. This paper 
examines the limits of the linear theory when the 
unsteady vortical flow velocity component is, in fact, 
not too small compared to the steady mean flow 
velocity. Thus, the interaction problem is investigated 
by imposing a finite-amplitude gust on the mean flow, 
and a nonlinear flow 


*Research Associate 
t Assistant Professor, Member of AXAA 
$ Professor, Associate Fellow of AIAA 
§Senior Research Scientist, Senior Member of AIAA 
^Senior Research Associate, Member of AIAA 
#Chief Scientist, Fellow of AIAA 

I 

American Institute of Aeroj 


analysis is then applied to examine possible deviations 
from the Imear theory predictions. The problem of gust- 
airfoil interaction is selected for this study as it presents 
a widely used benchmark in unsteady aerodynamics. 

To the authors’ best knowledge, no detailed, systematic 
study of high-intensity . gust response has been 
conducted so far (at least for realistic airfoil profiles), 
although nonlinear Euler and Navier-Stokes solvers 
have been in use for unsteady problems in 
turbomachinery for at least a decade (see [1] for a good 
review of the subject). A fully nonlinear viscous study 
of gust-airfoil interaction has been previously 
conducted by Lockard and Morris [2] with assumption 
of a small-amplitude perturbation. Hence, most 
observed nonlinear effects in [2] were related to gust- 
triggered instabilities in the viscous wake. In fact, in the 
present study, a similar mechanism giving rise to 
inviscid wake instability remains a crucial factor in 
describing the nonlinear response effects. In general, 
such effects to be investigated may include (e.g., [3], 
pp.863-873) generation of higher harmonics (both for 
single and multiple frequency input) and combination 
tones (for multiple-frequency input) resulting from 
nonlinear self-interaction and exchange of energy 
between multiple imposed frequencies. In certain cases, 
the effects of nonlinear interaction between the imposed 
harmonic(s) and the mean flow, including exchange of 
energy, can also be present. It is expected that these 
effects combined should become increasingly more 
noticeable at higher gust intensities when they would 
thus produce a more significant impact on the unsteady 
airfoil response, particularly on the unsteady 
aerodynamic loading and the acoustic directivity 
pattern. 

The present work extends the analysis of the gust-airfoil 
interaction problem first presented in [4] and [5]. 
Results of the nonlinear inviscid analysis based on 
explicit time marching with prefactored sixth-order 
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compact scheme for spatial differences, are compared 
with solutions from a frequency-domain linearized 
Euler solver GUST3D [6]. The unsteady response to 
one-dimensional and two-dimensional gusts is 
examined for both loaded and unloaded thick 
Joukowski airfoils in a subsonic flow. 

FORMULATION 


The analysis of nonlinear inviscid gust-airfoil 
interaction is based on the numerical solution to the 
nonlinear Euler equations which can be written in two- 
dimensional Cartesian coordinates as 

dQ dE dF n 

= 0 , ( 1 ) 


\Jl LA/U Kjy 


where Q is the solution vector, and F and G are the flux 
vectors: 
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and 


p = (y- 1) 
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Since the gust response is investigated for the complex- 
shaped Joukowksi airfoil, it requires recasting the 
equations in generalized curvilinear coordinates, with 
the chain-rule curvilinear Euler equations written as 


explicit boundary stencils for spatial derivatives [8]. A 
10* order explicit filter is used at every stage of the 
Runge-Kutta solver to provide dissipation. 

For the analysis of the gust response problem, a C-grid 
2D topology is used for the 12% thick Joukowski 
airfoil, with the grid extending at least 10 chords away 
from the surface in each direction. All the meshes are 
generated as single blocks by the GridPro commercial 
software [9]. More details on the numerical 
implementation are included below. 

GUST AND AIRFOIL CONFIGURATIONS AND 
FLOW INITIALIZATION 

In the current work, the nonlinear gust response is 
examined for a series of imposed gust intensities and 
frequencies. Each vortical gust harmonic is initially 
imposed on the mean flow with the following 
distribution: 


cos(cct + fy - at) 

(5) 

cos (ox + fiy — at) 


where e is the gust intensity relative to the mean flow, 
a and (3 are the gust wave numbers in the x and y 
directions, and co is the imposed gust frequency. The 
mean flow is defined far upstream from the airfoil as: 
p=l 

u - M 

v = 0 (6) 

P=~(P Y ) 

• • Y 

where M is the upstream mean flow Mach number, and 
y=1.4. For this paper’s computations, the Mach number 
is fixed at M- 0.5. 
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dQ d% dE dr] dE d% dF dr] dF 

— + — — + — - — + — + — — = 0 (4) 

dt dx d£ dx dr] dy dE, dy dr] 

The numerical formulation established in [4] and [5] 
uses the low-storage 4 th order 5-6 Low Dispersion and 
Dissipation Runge-Kutta scheme [7] for time marching, 
and the prefactored 6 th order compact scheme and 


The flow is initialized throughout the physical domain, 
with the gust superimposed on the mean flow. In 
addition, at the inflow 7 boundary 7 , the acoustic radiation 
condition given by Tam and Webb [10] is applied to the 
outgoing perturbation which is defined, e.g., for u- 

component as u BC = u boundary —u— U gust . 
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The following two configurations for the airfoil 
geometry are considered: 

(Al) Unloaded symmetric 12%-thick Joukowski airfoil 
(no angle of attack); 

(A2) Loaded 12%-thick 2%-cambered Joukowski 
airfoil at a two-degree angle of attack. 

Two configurations for the impinging gust are 
examined: 

(Gl) One-dimensional (ID) transverse gust with 

a = 2k,J3 =0,co = 2kM in (5); 

(G2) Two-dimensional (2D) gust with 
a = 2k, (5 = 2k, (O = 2kM in (5), 

where k is the reduced frequency of the gust non- 
dimensionalized by the half-chord and the upstream 
flow velocity, In the computations, both the gust 
intensity and reduced frequency are varied. 

RESULTS AND DISCUSSION 

The focus of this study is on nonlinear effects of high- 
intensity gust interaction with an airfoil. Hence, it is 
important to provide comparison for different gust 
intensity levels, with the corresponding low-intensity 
results assumed to fall within the linear response limits. 
At first, the results of the numerical analysis for all gust 
and airfoil configurations are validated against the 
GUST3D linear analysis calculations. Following this, a 
systematic study of gust intensity effects on various 
nonlinear response aspects is carried out, with the most 
important findings discussed and illustrated below. 

Numerical Implementation 

The computational code performing nonlinear unsteady 
analysis is written in Fortran 90 and uses the MPI 
library for parallelization. Numerical solutions are 
obtained through parallel implementation on the 
Beowulf-type Linux Athlon cluster consisting of eight 
dual-processor 1.4GElz AMD compute nodes," each with 
1GB of PC2100DDR memory. The nodes are 
connected by five full duplex fast ethemet connections 
to the other nodes in the cluster using a symmetric 
interconnection scheme. The MPICH implementation 
of MPI is used [11], 

The parallel performance of the code is shown in Figure 
1, with the noticeable decrease in the computational 
efficiency observed when more then eight processors 
are involved in the computations. This is due to the fact 
that the cluster starts to utilize both processors of the 
single nodes sharing the same bus allocation on the 
motherboard. For instance, the same computation 
performed with two processors from different nodes is 
observed to have the parallel efficiency of 95 %, 


compared to 73 % when two processors from the same 
node are used. This trend is typical of communication- 
intensive numerical implementations. 



Figure 1: Observed speedup for parallel runs. 

In this work, the simulations are run for longer time 
periods compared to [4] since stable, time-accurate 
solutions are sought not only for imposed gust 
frequencies but also for higher harmonics and 
combination tones. In terms of the non-dimensional 
time (normalized by the airfoil half-chord and the speed 
of sound), the values between 360 and 450 are found 
necessary for the more slowly converging loaded airfoil 
simulations to accurately resolve the mean flowfield 
and the acoustic field. More time iterations are required 
for lower gust frequency cases in order to march 
through an adequate number of longer periods. Some 
small differences in results compared to the similar 
cases in [4] can be attributed to longer convergence 
histories in the present runs. Finally, to quantify the 
computational effort, a 16-hour numerical run is 
required to march to the non-dimensional time value of 
360, which corresponds to 70,000 iterations on the 
600x240 grid using all sixteen processors. 

Validation 

The validation of the nonlinear code results with the 
linear solutions from the GUST3D solver has been 
previously accomplished in. [4], where thorough grid 
density and domain size sensitivity studies have also 
been conducted. The reader is referred to this paper for 
more details. In the present work, this validation is 
further extended to numerical simulations on meshes 
with grid densities 600x240 (for the cambered airfoil 
cases) and 591x238 (for the symmetric airfoil cases). 
Similar to [4], the grids are generated by GridPro using 
algebraic clustering around the profile in the normal 
direction (An=0.01) and near the trailing edge 


3 

American Institute of Aeronautics and Astronautics 



(Ax=0.01), with the stretching ratio of 1.05 used to 
expand to the far field. Note that some oscillations and 
loss of accuracy are found near the trailing edge, which 
is due to the singularity treatment in the grid topology. 
Nevertheless, a careful clustering of grid points near the 
trailing edge may help to keep these effects localized, 
with minimal impact on the overall solution. 

The results are presented for two reduced frequencies of 
the gust, k=0.1 in Figures 2 and 3, and k=1.0 in Figures 
4 and 5. For the selected airfoil and gust configurations 
(Al-Gl), (A1-G2), (A2-G1), (A2-G2), Figures 2 and 4 
show the acoustic intensity of the unsteady response at 
a distance of one chord length away from the 
centerpoint of the airfoil, while Figures 3 and 5 show 
results for the same quantity at a distance of four chord 
lengths. For each case, three results are compared 
including the low-intensity gust response (for e=0.02), 
the high-intensity gust response (for £=0.20), and the 
GUST3D linear solution [6,12], 

1 D Gust 

£= 0.02 




Figure 2:Disribution 


A detailed discussion of the high-intensity gust results 
is presented in the following paragraphs. The low- 
intensity solutions match well with GUST3D 
predictions, with only the case of high frequency and 
cambered profile showing apparent discrepancy (also in 
line with comparison in [4]). 

In general, the comparison is better for the unloaded 
airfoil cases. It should be also mentioned that the 
unsteady pressure distribution on the airfoil shows an 
excellent match for all the considered cases, with 
results nearly identical to the previous comparison in 
[4]. 


2D Gust 

£= 0.02 



,2 | atR=l (k=0.1) 
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1 D Gust 
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2D Gust 

e=0.02 




Figure 5:Disribution of |p’ 2 | at R=4 (k=1.0) 


Single-Frequency Gust 

Figures 2-5 also introduce a comparison between the 
solutions obtained for the low (£=0.02) and high (e=0.2) 
gust disturbance amplitudes. All the results in these and 
the following figures are scaled to the gust amplitude 
e=0.02 for proper comparison. Note that at the low 
reduced frequency of k=0.1, the directivity patterns for 
all the cases retain the same shape, but, at least in the 
near field for R=l, show a notable overall increase in 
the lobe acoustic intensity magnitudes, less evident in 
the far field (R=4). For the higher reduced frequency of 
k=1.0, the situation actually reverses, with only small 
lobe variation noticeable in the near field but rather 
significant change in the directivity shape in the far 
field, especially for the 2D gust cases. To examine this 
effect in more detail. Figure 6 shows the variation of 
the acoustic intensity in the far field over a range of 2D 
gust amplitudes for the loaded airfoil case. 


Note that multiple lobes start to develop in the 
directivity pattern with increasing gust amplitude, 
which may point to the appearance of new non-compact 
sources in the domain. Indeed, the numerical analysis 
indicates that this phenomenon is connected to the 
development of the inviscid wake instability triggered 
by the breakdown of the gust at the airfoil surface (first 
noticed by Lockard and Morris [2] for the viscous 
case). 

To illustrate this phenomenon. Figures 7-10 show 
contour plots for the instantaneous vertical velocity 
component. The results imply an increasing vorticity 
downstream in the wake, with the growth much more 
rapid for higher £. 
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Moreover, a similar trend is noted in Figure 11 for the 
RMS unsteady pressure fluctuations in the wake, which 
gives rise to the spike in the directivity pattern across 
the wake region, observed in Figure 6. 


harmonics is greatly increased at the higher gust 
intensity. A single case of unsteady response of the 
cambered airfoil to the ID gust with amplitude s=0.20 
also reveals harmonics of an unexplained spurious 
frequency. 



Figure 6: Effect of c on |p ,2 | Disribution at R=4 (k=1.0, 
2D gust, cambered airfoil) 
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Thus, the lobe distortion in Figure 6 is apparently 
attributed to the wave-like instability development in 
the wake as the radiation pattern from such instability 
waves may develop a multi-lobe response (this 
phenomenon is examined in [13] using an integral 
Kirchhoff technique). Also, from Figures 3 and 5, one 
may notice that the 2D gust produces a more significant 
nonlinear effect on directivity, which may possibly be 
attributed to the fact that the longitudinal gust 
component induces additional shear in the wake thus 
contributing to instability development. All these 
hypotheses will be validated in a future analysis. Note 
that the near-field results (R=l) do not reveal these 
directivity distortions since major no nlin ear wake 
evolution occurs further downstream. 

Figure 12 illustrates the contours of the calculated 
instantaneous unsteady pressure distribution in the 
domain. In order to examine the higher harmonics of 
the unsteady gust response, the FFT of the unsteady lift 
coefficient is carried out for all the test cases. Figure 13 
shows the results scaled to the same input gust 
amplitude for comparison. In the analysis, the FFT is 
performed over a period corresponding to the non- 
dimensional time of 407t. For this period, the code 
switches from an adjustable time step (satisfying the 
stability criteria) to the fixed one, using small steps to 
march through more than 30,000 iterations. It is 
obvious from Figure 13 that the level of higher 


Figure 7: Instantaneous velocity, Y-component: k=1.0, 
ID gust, e=0.02, symmetric airfoil. 



J L 



Figure 8: Instantaneous velocity, Y-component: k=1.0, 
ID gust, s=0.20, symmetric airfoil. 
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2D gust, s=0.02, cambered airfoil. x 

Figure 12: Contours of instantaneous unsteady pressure 
response. 2D gust, k=1.0, e=0.20, cambered airfoil. 


It should be mentioned that for gust intensities higher 
than s=0.2 and for higher gust frequencies, the 
nonlinear simulations start failing to converge as the 
spurious reflections start building up in the 
computational domain. This is attributed to the 
currently used formulation for the acoustic radiation 
boundary conditions based on the linear theory and the 
point source assumption. Future work will focus on 
resolving this issue to extend the analysis. This aspect is 
especially important because of the need to examine the 
nonlinear response at higher gust frequencies more 
typical of rotor/stator unsteady interactions. 

Figure 10: Instantaneous velocity, Y-component: k=1.0, 

2D gust, 8=0.20, cambered airfoil. Multiple-Frequency Gust 

In order to further investigate the nonlinear response at 
high gust amplitudes, several simulations have been run 
with superposition of two different gust frequencies. 
Only the loaded airfoil is considered. In Figure 14, the 

lOTIE-d ^ , 

FFT of the lift coefficient is presented for two 

j§|§ superimposed 2D gust frequencies with k=0.6 and 

k=1.0. The gusts are imposed with equal amplitudes, 
>£15 and comparison is presented for two cases with 8=0.06 
and 8=0.1. As expected, the unsteady response reveals 
the presence of combination tones with k=0.4 and 
k=1.6, in addition to higher harmonics. All of them 
x reveal higher amplitudes at higher gust intensities. 

Figure 1 1 : RMS unsteady pressure response contours Fmall y’ for ^ same lm P osed BPSt amplitudes, Figure 

downstream of airfoil (with R=4 line shown). 2D gust, 15 compares the unsteady responses due to the 

k=1.0, 8=0.20, cambered airfoil. superposition of ID vs. 2D gusts. Clearly, the 

superposition of 2D gusts generates higher levels of 
nonlinear response. 
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£= 0.02 


£= 0.20 




Figurel3: FFT of the lift coefficient for the single imposed gust frequency. 



Figure 14: FFT of the lift coefficient Figure 15: EFT of the lift coefficient for 

for k=0.6+1.0. Cambered airfoil, 2D k=0. 6+1.0. Cambered airfoil, ID and 2D 

gusts, cases with s=0.06 and 8=0.1 gust cases, 8=0.1 
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CONCLUSIONS 

This study focused on nonlinear effects of high- 
intensity gust interaction with an airfoil. To investigate 
the limits of validity for the Unearned approach, results 
of the nonlinear inviscid analysis were compared with 
the linearized Euler solutions in the parametric space of 
gust dimension, intensity, and frequency. Both loaded 
and unloaded airfoils in a subsonic flow were 
considered. 

Among the nonlinear effects, generation of higher 
harmonics and combination tones associated with 
nonlinear self-interaction and exchange of energy 
between multiple imposed frequencies was observed for 
gust ampHtudes approaching 20% of the mean flow. At 
low reduced frequency of k=0.1, the double-lobe 
directivity patterns tended to retain its shape but 
increased the overall acoustic intensity magnitudes. On 
the other hand, for the higher reduced frequency of 
k=1.0, most significant changes occurred in the 
directivity shapes starting to develop multiple lobes. 
r This effect was attributed to the development of 
inviscid wake instability acting as a new non-compact 
source. 

Overall, no significant deviation from the small- 
perturbation predictions was observed for gust 
amplitudes up to about 10% of the mean flow. It was 
further determined that, compared to the ID transverse 
gust, the 2D gust produced more significant nonlinear 
effects on the acoustic directivity in all cases 
considered, with the difference more pronounced at the 
higher reduced frequency. Nevertheless, no clear 
conclusion on the effects of gust frequency and airfoil 
loading on nonlinear response could be derived at this 
point. Future work is planned to extend this analysis to 
higher frequencies more typical of unsteady rotor-stator 
interactions, which will require development of a new 
approach to the acoustic boundary condition treatment 
more adequate for nonlinear computations. Effects of 
viscosity on nonlinear response particularly in the wake 
region will be examined with introduction of implicit 
time marching in the existing structure of the Very 
Large Eddy Simulation code used in this paper. The 
remaining discrepancies between the nonlinear analysis 
at low gust intensity and linear GUST3D solutions will 
be further investigated. 
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